library(tidyverse)
── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.4 ✔ tidyr 1.3.1
✔ purrr 1.0.2 ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(Seurat)
Loading required package: SeuratObject
Loading required package: sp
‘SeuratObject’ was built under R 4.4.0 but the current version is 4.4.1; it is recomended that you reinstall ‘SeuratObject’ as the ABI for R may have changed
‘SeuratObject’ was built with package ‘Matrix’ 1.7.0 but the current version is 1.7.1; it is recomended that you reinstall ‘SeuratObject’ as the ABI for ‘Matrix’ may have changed
Attaching package: ‘SeuratObject’
The following objects are masked from ‘package:base’:
intersect, t
library(patchwork)
library(SingleR)
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘matrixStats’
The following object is masked from ‘package:dplyr’:
count
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps,
colMadDiffs, colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates,
colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians, colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse,
rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps, rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads, rowWeightedMeans,
rowWeightedMedians, rowWeightedSds, rowWeightedVars
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Attaching package: ‘BiocGenerics’
The following object is masked from ‘package:SeuratObject’:
intersect
The following objects are masked from ‘package:lubridate’:
intersect, setdiff, union
The following objects are masked from ‘package:dplyr’:
combine, intersect, setdiff, union
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, intersect, is.unsorted,
lapply, Map, mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, setdiff, table, tapply, union, unique, unsplit,
which.max, which.min
Loading required package: S4Vectors
Attaching package: ‘S4Vectors’
The following objects are masked from ‘package:lubridate’:
second, second<-
The following objects are masked from ‘package:dplyr’:
first, rename
The following object is masked from ‘package:tidyr’:
expand
The following object is masked from ‘package:utils’:
findMatches
The following objects are masked from ‘package:base’:
expand.grid, I, unname
Loading required package: IRanges
Attaching package: ‘IRanges’
The following object is masked from ‘package:sp’:
%over%
The following object is masked from ‘package:lubridate’:
%within%
The following objects are masked from ‘package:dplyr’:
collapse, desc, slice
The following object is masked from ‘package:purrr’:
reduce
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
Attaching package: ‘SummarizedExperiment’
The following object is masked from ‘package:Seurat’:
Assays
The following object is masked from ‘package:SeuratObject’:
Assays
library(celldex)
Attaching package: ‘celldex’
The following objects are masked from ‘package:SingleR’:
BlueprintEncodeData, DatabaseImmuneCellExpressionData, HumanPrimaryCellAtlasData, ImmGenData, MonacoImmuneData, MouseRNAseqData, NovershternHematopoieticData
library(pheatmap)
# Speeds up scTransform
#BiocManager::install("glmGamPoi")
# Load AR and Tol datasets
AR_data <- CreateSeuratObject(Read10X(data.dir = "~/Desktop/R_projects/Dr_Krummey_Lab/heartTransplant_data/GSE262851_RAW/AR/"), min.cells = 3, min.features = 200)
AR_data$condition <- "AR"
AR_data <- PercentageFeatureSet(AR_data, pattern = "^mt-", col.name = "percent.mt")
VlnPlot(AR_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
# I filtered a bit looser, nFeatureRNA > 5000 vs 2500, percent.mt < 15 vs 5
AR_data <- subset(AR_data, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
VlnPlot(AR_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
Tol_data <- CreateSeuratObject(Read10X(data.dir = "~/Desktop/R_projects/Dr_Krummey_Lab/heartTransplant_data/GSE262851_RAW/Tol/"), min.cells = 3, min.features = 200)
Tol_data$condition <- "Tol"
Tol_data <- PercentageFeatureSet(Tol_data, pattern = "^mt-", col.name = "percent.mt")
VlnPlot(Tol_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
Tol_data <- subset(Tol_data, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)
VlnPlot(Tol_data, features = c("nCount_RNA", "nFeature_RNA", "percent.mt"))
Warning: Default search for "data" layer in "RNA" assay yielded no results; utilizing "counts" layer instead.
# Lets merge asap so we know the same transformations are applied
# Project is just a label to group the merged datasets and I never reference it. Feel free to change!
data.merge <- merge(AR_data, Tol_data, add.cell.ids = c("AR", "Tol"), project = "Heart")
# SCTransform instead of normalize -> variable features -> scale. More up-to-date
data.merge <- SCTransform(data.merge, vars.to.regress = "percent.mt", verbose = FALSE)
data.merge <- RunPCA(data.merge, verbose = FALSE)
# dims of PCA to choose. Probably some fine tuning here for dims and number of neighbors, but I think it looks fine enough already
# It is always good to be thinking whether you are overfitting or underfitting data though!
data.merge = FindNeighbors(data.merge, dims = 1:30, reduction = "pca")
Computing nearest neighbor graph
Computing SNN
data.merge = FindClusters(data.merge, resolution = 0.5, cluster.name = "unintegrated_clusters")
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5321
Number of edges: 181553
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8922
Number of communities: 17
Elapsed time: 0 seconds
data.merge <- RunUMAP(data.merge, dims = 1:30, reduction = "pca", reduction.name = "umap.unintegrated", verbose = FALSE)
Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session
# Always visualize!
# The plot is showing batch effects, so we definitely need to integrate
DimPlot(data.merge, reduction = "umap.unintegrated", group.by = c("condition", "unintegrated_clusters"), combine = FALSE)
[[1]]
[[2]]
# Now we integrate the two layers using the Harmony algorithm with the PCA we calculated
data.integrated <- IntegrateLayers(object = data.merge, method = HarmonyIntegration, orig.reduction = "pca", new.reduction = "harmony", normalization.method = "SCT", verbose = FALSE)
Warning: HarmonyMatrix is deprecated and will be removed in the future from the API in the futureWarning: Warning: The parameters do_pca and npcs are deprecated. They will be ignored for this function call and please remove parameters do_pca and npcs and pass to harmony cell_embeddings directly.Warning: Warning: The parameter tau is deprecated. It will be ignored for this function call and please remove parameter tau in future function calls. Advanced users can set value of parameter tau by using parameter .options and function harmony_options().Warning: Warning: The parameter block.size is deprecated. It will be ignored for this function call and please remove parameter block.size in future function calls. Advanced users can set value of parameter block.size by using parameter .options and function harmony_options().Warning: Warning: The parameter max.iter.harmony is replaced with parameter max_iter. It will be ignored for this function call and please use parameter max_iter in future function calls.Warning: Warning: The parameter max.iter.cluster is deprecated. It will be ignored for this function call and please remove parameter max.iter.cluster in future function calls. Advanced users can set value of parameter max.iter.cluster by using parameter .options and function harmony_options().Warning: Warning: The parameter epsilon.cluster is deprecated. It will be ignored for this function call and please remove parameter epsilon.cluster in future function calls. Advanced users can set value of parameter epsilon.cluster by using parameter .options and function harmony_options().Warning: Warning: The parameter epsilon.harmony is deprecated. It will be ignored for this function call and please remove parameter epsilon.harmony in future function calls. If users want to control if harmony would stop early or not, use parameter early_stop. Advanced users can set value of parameter epsilon.harmony by using parameter .options and function harmony_options().
# After integrating we re-find our neighbors and clusters and run UMAP to check our batch effects visually
data.integrated <- FindNeighbors(data.integrated, reduction = "harmony", dims = 1:30)
Computing nearest neighbor graph
Computing SNN
data.integrated <- FindClusters(data.integrated, resolution = 0.5, cluster.name = "harmony_clusters")
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 5321
Number of edges: 188646
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8902
Number of communities: 14
Elapsed time: 0 seconds
data.integrated <- RunUMAP(data.integrated, reduction = "harmony", dims = 1:30, reduction.name = "umap.harmony", verbose = FALSE)
# Much nicer, AR and TOL form well-mixed clusters
DimPlot(data.integrated, reduction = "umap.harmony", group.by = c("condition", "harmony_clusters"), combine = FALSE)
[[1]]
[[2]]
# I would suggest looking into Azimuth instead of SingleR. Better fine labeling
ref <- celldex::ImmGenData()
data.integrated_fine_pred <- SingleR(test = GetAssayData(data.integrated, layer = 'counts'),
ref = ref,
labels = ref$label.fine)
data.integrated_main_pred <- SingleR(test = GetAssayData(data.integrated, layer = 'counts'),
ref = ref,
labels = ref$label.main)
data.integrated$singleR.labels.fine <- data.integrated_fine_pred$labels[match(rownames(data.integrated@meta.data), rownames(data.integrated_fine_pred))]
data.integrated$singleR.labels.main <- data.integrated_main_pred$labels[match(rownames(data.integrated@meta.data), rownames(data.integrated_main_pred))]
DimPlot(data.integrated, reduction = "umap.harmony", group.by = "singleR.labels.main")
DimPlot(data.integrated, reduction = "umap.harmony", group.by = "singleR.labels.fine")
head(sort(table(data.integrated_fine_pred$labels), decreasing=TRUE))
Macrophages (MF.11CLOSER.SALM3) Monocytes (MO.6C+II-) T cells (T.8EFF.OT1.D10LIS) Macrophages (MF.103-11B+24-) Endothelial cells (BEC) Macrophages (MF.II+480LO)
1408 561 298 228 219 204
head(sort(table(data.integrated_main_pred$labels), decreasing=TRUE))
Macrophages T cells Monocytes NKT B cells Endothelial cells
2807 826 352 334 271 220
# 826 T cells
# Finding more ways to plot our results, checking annotation of cells and pruned cells
plotScoreHeatmap(data.integrated_main_pred)
plotDeltaDistribution(data.integrated_main_pred)
pattern <- "T cells \\(T\\.(8|CD8).*\\)"
# Shorter way to add metadata
data.integrated$is_cd8 <- grepl(pattern, data.integrated@meta.data$singleR.labels.fine)
DimPlot(data.integrated, reduction = "umap.harmony", group.by = "is_cd8")
cd8_subset <- subset(data.integrated, subset = is_cd8 == TRUE)
# Rerun SCTransform on the subset, since the global structure of variance changed
cd8_subset <- SCTransform(cd8_subset, vars.to.regress = "percent.mt", verbose = FALSE)
Warning: Different cells and/or features from existing assay SCT
# Sub-clustering of CD8 cells
cd8_subset <- RunPCA(cd8_subset, reduction.name = "cd8.pca", verbose = FALSE)
Warning: Key ‘PC_’ taken, using ‘cd8pca_’ instead
cd8_subset <- FindNeighbors(cd8_subset, reduction = "cd8.pca", dims = 1:30)
Computing nearest neighbor graph
Computing SNN
cd8_subset <- FindClusters(cd8_subset, resolution = 0.5, cluster.name = "cd8_clusters")
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
Number of nodes: 1061
Number of edges: 43975
Running Louvain algorithm...
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.6805
Number of communities: 5
Elapsed time: 0 seconds
cd8_subset <- RunUMAP(cd8_subset, dims = 1:30, reduction = "cd8.pca", verbose = FALSE, reduction.name = "cd8.umap")
DimPlot(cd8_subset, reduction = "cd8.umap", group.by = c("condition", "cd8_clusters"), combine = FALSE)
[[1]]
[[2]]
# Finding DE genes
Idents(cd8_subset) <- "cd8_clusters"
cd8_subset <- PrepSCTFindMarkers(cd8_subset)
Found 2 SCT models. Recorrecting SCT counts using minimum median counts: 5574
| | 0 % ~calculating
|+++++++++++++++++++++++++ | 50% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s
tcell.cd8.markers <- FindMarkers(cd8_subset, ident.1 = 0, ident.2 = 2)
For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the presto package
--------------------------------------------
install.packages('devtools')
devtools::install_github('immunogenomics/presto')
--------------------------------------------
After installation of presto, Seurat will automatically use the more
efficient implementation (no further action necessary).
This message will be shown once per session
| | 0 % ~calculating
|+ | 1 % ~15s
|+ | 2 % ~14s
|++ | 3 % ~14s
|++ | 4 % ~14s
|+++ | 5 % ~14s
|+++ | 6 % ~14s
|++++ | 7 % ~14s
|++++ | 8 % ~13s
|+++++ | 9 % ~13s
|+++++ | 10% ~13s
|++++++ | 11% ~13s
|++++++ | 12% ~13s
|+++++++ | 13% ~13s
|+++++++ | 14% ~13s
|++++++++ | 15% ~13s
|++++++++ | 16% ~13s
|+++++++++ | 17% ~12s
|+++++++++ | 18% ~12s
|++++++++++ | 19% ~12s
|++++++++++ | 20% ~12s
|+++++++++++ | 21% ~12s
|+++++++++++ | 22% ~12s
|++++++++++++ | 23% ~12s
|++++++++++++ | 24% ~11s
|+++++++++++++ | 25% ~11s
|+++++++++++++ | 26% ~11s
|++++++++++++++ | 27% ~11s
|++++++++++++++ | 28% ~11s
|+++++++++++++++ | 29% ~11s
|+++++++++++++++ | 30% ~11s
|++++++++++++++++ | 31% ~10s
|++++++++++++++++ | 32% ~10s
|+++++++++++++++++ | 33% ~10s
|+++++++++++++++++ | 34% ~10s
|++++++++++++++++++ | 35% ~10s
|++++++++++++++++++ | 36% ~10s
|+++++++++++++++++++ | 37% ~09s
|+++++++++++++++++++ | 38% ~09s
|++++++++++++++++++++ | 39% ~09s
|++++++++++++++++++++ | 40% ~09s
|+++++++++++++++++++++ | 41% ~09s
|+++++++++++++++++++++ | 42% ~09s
|++++++++++++++++++++++ | 43% ~09s
|++++++++++++++++++++++ | 44% ~08s
|+++++++++++++++++++++++ | 45% ~08s
|+++++++++++++++++++++++ | 46% ~08s
|++++++++++++++++++++++++ | 47% ~08s
|++++++++++++++++++++++++ | 48% ~08s
|+++++++++++++++++++++++++ | 49% ~08s
|+++++++++++++++++++++++++ | 50% ~07s
|++++++++++++++++++++++++++ | 51% ~07s
|++++++++++++++++++++++++++ | 52% ~07s
|+++++++++++++++++++++++++++ | 53% ~07s
|+++++++++++++++++++++++++++ | 54% ~07s
|++++++++++++++++++++++++++++ | 55% ~07s
|++++++++++++++++++++++++++++ | 56% ~07s
|+++++++++++++++++++++++++++++ | 57% ~06s
|+++++++++++++++++++++++++++++ | 58% ~06s
|++++++++++++++++++++++++++++++ | 59% ~06s
|++++++++++++++++++++++++++++++ | 60% ~06s
|+++++++++++++++++++++++++++++++ | 61% ~06s
|+++++++++++++++++++++++++++++++ | 62% ~06s
|++++++++++++++++++++++++++++++++ | 63% ~06s
|++++++++++++++++++++++++++++++++ | 64% ~05s
|+++++++++++++++++++++++++++++++++ | 65% ~05s
|+++++++++++++++++++++++++++++++++ | 66% ~05s
|++++++++++++++++++++++++++++++++++ | 67% ~05s
|++++++++++++++++++++++++++++++++++ | 68% ~05s
|+++++++++++++++++++++++++++++++++++ | 69% ~05s
|+++++++++++++++++++++++++++++++++++ | 70% ~04s
|++++++++++++++++++++++++++++++++++++ | 71% ~04s
|++++++++++++++++++++++++++++++++++++ | 72% ~04s
|+++++++++++++++++++++++++++++++++++++ | 73% ~04s
|+++++++++++++++++++++++++++++++++++++ | 74% ~04s
|++++++++++++++++++++++++++++++++++++++ | 75% ~04s
|++++++++++++++++++++++++++++++++++++++ | 76% ~04s
|+++++++++++++++++++++++++++++++++++++++ | 77% ~03s
|+++++++++++++++++++++++++++++++++++++++ | 78% ~03s
|++++++++++++++++++++++++++++++++++++++++ | 79% ~03s
|++++++++++++++++++++++++++++++++++++++++ | 80% ~03s
|+++++++++++++++++++++++++++++++++++++++++ | 81% ~03s
|+++++++++++++++++++++++++++++++++++++++++ | 82% ~03s
|++++++++++++++++++++++++++++++++++++++++++ | 83% ~03s
|++++++++++++++++++++++++++++++++++++++++++ | 84% ~02s
|+++++++++++++++++++++++++++++++++++++++++++ | 85% ~02s
|+++++++++++++++++++++++++++++++++++++++++++ | 86% ~02s
|++++++++++++++++++++++++++++++++++++++++++++ | 87% ~02s
|++++++++++++++++++++++++++++++++++++++++++++ | 88% ~02s
|+++++++++++++++++++++++++++++++++++++++++++++ | 89% ~02s
|+++++++++++++++++++++++++++++++++++++++++++++ | 90% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++ | 91% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++ | 92% ~01s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 93% ~01s
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 95% ~01s
|++++++++++++++++++++++++++++++++++++++++++++++++ | 96% ~01s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97% ~00s
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 98% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 99% ~00s
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=15s
tcell.cd8.markers
FeaturePlot(cd8_subset, features = "Gcnt1", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Gzmb", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Spn", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Sell", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Prf1", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Spn", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Cxcr3", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Havcr2", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Ifng", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Icos", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Prdm1", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Tox", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Gzmk", reduction = "cd8.umap")
#FeaturePlot(cd8_subset, features = "Gbma", reduction = "cd8.umap") NOT FOUND
#FeaturePlot(cd8_subset, features = "Gzmh", reduction = "cd8.umap") NOT FOUND
FeaturePlot(cd8_subset, features = "Gzmb", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Lag3", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Tigit", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Irf4", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Tcf7", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Pdcd1", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Tnf", reduction = "cd8.umap")
FeaturePlot(cd8_subset, features = "Tbx21", reduction = "cd8.umap")
#Save the plots as PNGs
library(ggplot2)
genes <- c(“Gcnt1”, “Gzmb”, “Spn”, “Sell”, “Prf1”, “Cxcr3”, “Havcr2”, “Ifng”, “Icos”, “Prdm1”, “Tox”, “Gzmk”, “Lag3”, “Tigit”, “Irf4”, “Tcf7”, “Pdcd1”, “Tnf”, “Tbx21”)
for (gene in genes) { plot <- FeaturePlot(cd8_subset, features = gene, reduction = “cd8.umap”) ggsave(filename = paste0(“CD8_”, gene, “.png”), plot = plot, width = 6, height = 4, dpi = 300) }
# From the list Dr. Krummey sent, find the differential expression of the genes of interest (done)
# Look back to just NON-CD4 T cells to see if there is a difference in that and solely using/looking at CD8+ T cells